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! ! Abstract 

Q-i The present paper suggests a method for obtaining incompressible solenoidal ve- 

locity vectors that satisfy approximately the desired immersed velocity boundary 
conditions. The method employs merely the mutual kinematic relations between 
the velocity and vorticity fields (i.e, the curl and Laplacian operators). An initial 
non-solenoidal velocity field is extended to a regular domain via a zero- velocity mar- 
gin, where an extended vorticity is found. Re-calculation of the velocities (subjected 
to appropriate boundary conditions), yields the desired solenoidal velocity vector. 
The method is applied to the two- and three-dimensional problems for the homo- 
geneous Dirichlet, as well as periodic boundary conditions. The results show that 
1-H the solenoidality is satisfied up to the machine accuracy for the periodic boundary 

■^j- conditions (employing the Fourier-spectral solution method), while an improvement 

in the solenoidality is achievable for the homogenous boundary conditions. 

• • 

J> Key words: solenoidal velocity vectors; immersed boundary conditions; 

velocity-vorticity relation; harmonic functions; periodic and homogenous 
S_i boundary conditions 



1 Introduction 



The Navier-Stokes Equations (NSE) are derived for a continuum. Therefore, 
many difficulties arise in dealing with any kind of discontinuity, including 
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finite-size solution domains (simply-connected or multiply-connected). Partic- 
ularly in an incompressible flow, the problem usually appears as non-solenoidal 
velocity fields. Consequently, satisfying the mass conservation has been one 
of the most challenging issues of almost all algorithms of incompressible flow 
simulations, whether in the primitive variables formulations or the other for- 
malisms. As it is revealed, the boundary conditions have a vital role here; and 
therefore, appropriate definition of the boundary conditions received signifi- 
cant attention in the literature [9J12|I13] . 

The literature of the subject is too rich to completely be summarized in this 
paper (see e.g. [13] and the references in there); however, as the most pertinent 
works, we shall make reference to the works of Quartapelle and Valz-Gris [UJ , 
and also Quartapelle [12], which studied the issue extensively. Calhun [3] used 
the immersed boundary technique in the vorticity-stream function formula- 
tion, and introduced an appropriate distribution of the vorticity in order to 
satisfy the overall mass balance. From another viewpoint, the problem was 
studied in relation to the compatibility conditions (also called the space-time 
corner singularity). The required compatibility conditions for the NSE was 
studied by Temam [TT], Gresho and Sani [PJ, and also by Boyd and Flyer [TJ, 
and then by Flyer and Swarztrauber [5]). 

The main contribution of the present paper is to suggest a method for con- 
struction of a nearly solenoidal velocity field, which satisfies approximately the 
desired immersed velocity boundary conditions, using the kinematic velocity- 
vorticity relation. The central idea is to use the well-known property that the 
divergence of the velocities, obtained from the velocity-vorticity equation is 
a harmonic function; and therefore, gets its maximum values on the bound- 
aries. Consequently, the divergence can be reduced (or even vanished) all over 
the field, by restricting (or vanishing) the divergence near the boundaries. In 
the present method, these restrictions are implemented by considering a zero- 
velocity margin in the vicinity of the outer regular boundaries. 
Embedding non-regular solution domains in a bigger regular domain was 
shown to be useful in solution of the scalar Poisson's problems [2lfH] . Now the 
present work may be seen as an extension of the method to the vector Pois- 
son's problems that the solenoidal solutions are desired. As it will be shown, 
a solenoidal velocity field (within the machine accuracy) is achievable for the 
periodic boundary conditions (that is, the no-boundary problems), and also 
improvements in the solenoidality is possible for the homogeneous Dirichlet 
boundary conditions. The method can be used in recently popularized one- 
fluid [19], or one-continuum [16] immersed boundary methods, as it was used 
successfully our previous work [To] . 

The paper is continued by description of the suggested algorithm, follows by 
a section contains our arguments about the reasons that the method works 
(that is, §[3]). As our numerical experiments, two- and three-dimensional fixed 
and moving boundary problems are examined, with homogeneous as well as 
periodic boundary conditions. 



2 



Zero velocity margin (£>) 



Fig. 1. The flow domain Oy, together with fixed or moving obstacle(s) O s , and 
other given- velocity regions f2 u , are embedded in the regular solution domain D 
(with regular boundary Tp), via a zero- velocity margin B. 



2 The suggested procedure 



Given an arbitrary (may be non-solenoidal) velocity vector u* = (w*, ■ ■ ■ , u^) G 
C 2 (Qf), denned on the flow domain Qf 6 M. d , with arbitrary boundary Tg, 
where d = 2, 3; one can find the corresponding vorticity vector u>* = V x u* 
on the Qf. Now, it is well-known that solution of the Poisson equations 



V 2 u = -V x u* 



(1) 



on Qf implies 

V 2 (V ■ u) = 0; 



(2) 



which means although V • u is not necessarily zero, it is a harmonic func- 
tion, and therefore, gets its maximum values on the boundaries. The above 
statement has been demonstrated several times by many authors [111112] . and 
motivated many attempts in obtaining solenoidal velocities by vanishing the 
divergence right at the boundaries 

Now, we suggest limiting the divergence inside the solution domain via con- 
trolling the divergence in the vicinity of the outer boundaries (not merely at 
the boundaries), via definition of a zero- velocity margin. According to Fig. [IJ 
The regular solution domain D = D + contains all the fluid-solid system, 
including the moving regions of the fluid with given velocity boundary condi- 
tions (that is fl u ), in addition to the fixed or moving immersed solid objects. 



Initial velocity u*(£2/) 
(may be non-solenoidal) 

1 

Extension to the regular domain 

u tD (D) 



f 

Vorticity on the extended domain 

u D = V x u* D 

1 

Solution of the Poisson's Eqns. 

V 2 U° = -V X U! D 

I 

Extraction of the desired velocities 

U (n / ) 



Fig. 2. Flowchart of the proposed algorithm. The initial velocities on the flow domain 
£lf is extended to a bigger regular domain D, where the extended vorticity co D is 
obtained and the Poisson's equations are solved on. Then the modified velocities u 
are extracted from u D . 

Now this fluid-solid system is embedded in a regular domain D via a zero- 
velocity margin B. 

With the above definitions in mind, the following procedure is suggested: 

i) Given the arbitrary (may be non-solenoidal) velocity vector u*(f2j) denned 
on Qf, the extended velocity u* D (D), defined on the regular domain D, 
is found via addition of the zero-velocity margin B, and other immersed 
regions (Q u and others) to the u*. 
ii) The vorticity u D = V x u* D is obtained on D. 

iii) The modified extended velocity vector u D is obtained from solution of the 
Poisson's equations 

V 2 u D = -Vxw° xeD (3) 

subjected to some boundary conditions on To, which will be discussed in 
the next section. 

iv) The desired modified velocity vector u(Q/) is obtained by discarding u D (B) 
and the other immersed regions (fl u and others) . 
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The above procedure is illustrated in Fig. |2j 

Remark 1 In fact for the multiply connected domains, the extremum of the 
divergence may occur on the boundaries of the holes (e.g. Q u or others). How- 
ever, in the present work we followed an immersed boundary approach, that 
is, u* D includes the velocities of all the holes of D. In this way, the extremum 
of the divergence occurs on Tp, and therefore, it is just needed to control the 
divergence on Tp. 



3 On the Solenoidality of u 



We shall examine the solenoidality of u. To this end, we consider two different 
boundary conditions for Eq. (|3]), that is, the periodic boundary conditions, 
which is suitable for the Fourier spectral solvers, and the homogenous Dirichlet 
boundary conditions, which can be implemented easily in the finite-difference 
or finite-volume solvers. 



3.1 Periodic Boundary Conditions 



In the case of periodic boundary conditions (that is, no-boundary problems), 
equation ^ can be transformed into the Fourier space, which leads to 

u D = ~i^r 9 x CP, (4) 



in which (•) = FT(-); and k = (ki, • • • , kd) stands for the wavenumber vector, 
and % = y/—l. Obviously, u D is perpendicular to the wavenumbers vector k, 
that is 

k ■ u D = 0, (5) 



which is a translation of solenoidality of u D in the Fourier space [3]. Now the 
desired (solenoidal) velocity vector u can be obtained simply by ignoring the 
solution in the B and the other immersed regions. More than the above ar- 
guments, our numerical experiments have also confirmed that in this case u 
is solenoidal within the machine accuracy. In fact, achieving perfect conserva- 
tion of mass for the no-boundary problems is not so surprising, because many 
assumptions that the NSE are based on are satisfied [7]. 
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3.2 Homogenous Dirichlet boundary conditions 



For the homogenous boundary conditions on Tp, the modified extended ve- 
locity u D is not strictly solenoidal. However, since V 2 (V • u D ) = 0, the diver- 
gence of u D gets its extremum on the boundary T^. In the other words, the 
extremum of divergence occurs in £>, and therefore, discarding the solution in 
B reduces the divergence such that 

HV-ulU < HV-u^Hoo. (6) 



In fact, the maximum values of V • u are discarded by cutting the solution 
in B\ and therefore, the divergence of u is less than the divergence of u D . 
As it will be seen in our numerical experiments, finite-difference solution of 
Eqns. ^ with homogeneous boundary conditions, modifies the solenoidality 
of an initially non-solenoidal velocity field. 



6 



1,25 




-1.25 1.25 

Fig. 3. Left: The modified Chandrashkar-Reid flow, defined on Clf = HxH^ < 1. 
The no-slip condition is satisfied formally on the boundaries Tg = {x, HxH^ = 1}. 
Right: The solution domain is extended to D = {x, HxH^ < 1.25}, via addition of 
the zero velocity margin B = D \ O. 



4 Numerical Experiments 

In this section, the proposed algorithm is assessed via examination of some 
fixed, as well as moving boundary problems in two and three dimensions. We 
implemented the method both in the Fourier-spectral, and the finite-difference 
solutions. 



4-1 The modified Chandrashkar-Reid flow 

In the study of the orthogonal basis functions that satisfy the no-slip con- 
ditions, Chandrashkar and Reid proposed a velocity field that satisfies the 
no-slip boundary condition in one direction, and is periodic in the other direc- 
tion [3]. Then this flow was appeared in the work of Rempfer |I3]. In the light 
of their works, it is convenient here to construct a velocity field, which satisfies 
the no-slip conditions in both directions. Tho this end, we first introduce a 
one-dimensional characteristic function 

7T 

<t>{xi) = cos(Axj) + A\ cosh(-Xi); x { e [-1, +1], (7) 



where A = 2.64, and A\ = — cos A/ cosh(7r/2) are chosen such that 0(±1) = 0, 
within the machine accuracy. Using this characteristic function, we construct 
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Fig. 4. Distribution of the divergence for the initial modified Chandrashkar-Reid 
flow. Although the velocities are obtained from a streamfunction, the divergence is 
not zero and its maximum values are occurred at the boundaries. 



the desired streamfunction 



in which the normalization factor G is chosen such that the unit maximum 
velocity occurs in the field, that is 

G=(l + A A )[^ A sinh(^)-Asin(^)]. (9) 



Now the velocities are directly obtainable from this streamfunction as 



d Id 

Ui = Q^-ip(xi, x 2 ) = -^<p{xi)-^<p{x 2 ), (10) 

d Id 



The left panel of Fig. [3] shows the above velocity filed on a 256 2 -point grid, 
which formally satisfies the no-slip conditions in both directions. In this con- 
figuration we defied Qf = \\x-\\oq < 1. 

It is worth mentioning that although the above velocity vector is obtained from 
a streamfunction, it is not solenoidal. In fact, the divergence can be written 
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as 



v..-I 



(12) 



which obviously is not necessarily zero on the boundaries. For example, on the 
left boundary <fi(x = — 1) = 0, while §^4>{x = — 1) ^ 0. 

This fundamental problem is one of the difficulties that the ip — u formulation 
is suffered from. In fact, in solution of V 2 ip = —w, two sets of boundary con- 
ditions are needed to be satisfied (that is, if)(Tf) = ipn and J^ip — 0), which 
makes the problem over-determined |13j . The divergence of the flow filed is 
shown in Figure. |4| As one can see, the divergences are not zero, especially 
right at the boundaries. 

Now the method is applied to this flow. To this end, the solution domain is 
extended to a bigger regular domain D = {HxH^ < 1.25}, and then the cor- 
responding vorticity u* is found on D. The extended velocity field is shown 
in the right panel of Fig. [3] Now Eqns. (|3]) are solved for the periodic as 
well as homogenous boundary conditions using the Fourier spectral, and the 
finite-difference methods; and the desired velocity u is extracted from u*. The 
distributions of V • u for both cases are shown in Fig. [5] As one can see, for 
the Fourier-spectral solutions the divergences are reduced to the machine ac- 
curacies (see the left panne of Fig. [5]). On the other hand, for the homogenous 
boundary conditions (which Eqns. ^ are solved using the finite-difference 
method), although the divergences are no longer vanished completely, how- 
ever, reduced appreciably. Table ?? summarizes the maximum values of the 
divergences in the above cases. 

It is worth mentioning that although the above procedure improves the 
solenoidality of the flow, it may change the velocity boundary conditions. 
In fact, the velocities are changed all over the field, especially near the bound- 
aries. This issue will be discussed in detail in the next section, because for the 
present test case boundary condition changes have not been sensible. 



4-2 Moving square cylinder 



In this section, the capability of the method in dealing with multiply-connected 
domains is examined. In this context, the practically important (and academ- 
ically challenging) problem of a moving obstacle in a quiescent fluid is con- 
sidered. Since the objective of the present paper is not implementation of the 
non-regular boundaries, a square cylinder, coinciding with the Cartesian grid, 
is considered. For the problem of non-regular immersed boundaries one can 
see the related references (see e.g. [HUS], or [2"Tj). 



9 



Fig. 5. The divergences for the extended modified Chandrashkar-Reid flow. Left: 
Fourier spectral solution for the periodic boundary conditions on Tj). The diver- 
gence is reduced to the machine accuracy. Right: Finite-difference solution for the 
homogenous boundary conditions on To (the dashed lines show the divergence of 
the extended velocity u D , and the solid lines show the desired final velocity u). 
The divergences are no longer within the machine accuracy, however, appreciable 
reduction is observable (c.f. Fig. El). 



u.-1 

□ 




Fig. 6. Left: Moving square cylinder in a quiescent fluid. The immersed velocity 
boundary conditions (u\, 112) = (0, 1) are implemented sharply on the square cylin- 
der Q, s = {\\xi — 7r||oo < 10/T287r}, placed in the flow domain [0, 2tt]. Right: The 
resulting divergences. As one can see, the continuity is lost because of implementa- 
tion of the boundary conditions. 

Consider a regular fluid-solid domain D, consists of a square cylinder Cl s = 
{Halloo — 10/1287r}, and the flow domain Qf = D — Cl s . Moreover, assume 
that it is desired to implement the velocity boundary condition u s = (1,0) 
to the square cylinder (which will be treated here as an immersed boundary 
condition). Implementation of such a boundary condition imposes some dis- 
continuities to the solution for t — > 0, the problem which is known as the 
corner singularity, or non-consistency of the boundary and initial conditions 
and were studied previously (e.g. see [lTlllTT8"f8"] ). The initial velocity and the 
distribution of the divergences are shown in Fig. k3} calculated on a 256 2 - point 
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Fig. 7. The modified velocities for the moving square cylinder with homogenous 
boundary conditions. Left: The induced velocities are observable almost on the 
whole of the regular domain D. Right: The maximum values of the divergences are 
reduced appreciably (c.f. Fig. [6]), and occur on the outer boundaries. 

grid. As one can see, the continuity is lost, especially on the boundaries of the 
moving obstacle, because of implementation of the sharp immersed velocity 
boundary conditions. 

Now the method is applied to this flow for the homogeneous and periodic 
boundary conditions. It should be noted that since the velocities are already 
zero in the vicinity of the regular boundary Tp, the zero velocity margin B is 
not needed indeed. First we consider the homogeneous boundary conditions. 
In this case, the vorticity u D is found, and the Poisson's equations ([3j are dis- 
cretized on D using the finite-difference method. The homogeneous Dirichlet 
boundary conditions are implemented to the discretized equations, and the 
resulting linear problem is solved using the point successive over relaxation 
(PSOR) method. The modified velocity field and the divergences are shown 
in Fig. [7| As one can see in the left panel, the initial immersed non-solenoidal 
velocities are now spread across the domain D, such that the non-zero veloci- 
ties are observable near the homogeneous outer boundaries (and therefore, the 
divergences are reduced near the square). As a consequence, the divergences 
in the right panel show an appreciable reduction in the maximum values (with 
the factor of ~ 1CF 4 in comparison to Fig. [6]), and furthermore, these maxi- 
mum values are now occurred on the outer (regular) boundaries. 
For the periodic boundary conditions, since the initial velocities u* are van- 
ishing in the vicinity of the regular boundary r^, the Fourier spectral solver 
may be used easily [2|14||15j . The vorticity (on D) is obtained in the Fourier 
space, and the modified velocities are found via solution of Eqns. The 
modified velocities and the divergences are shown in Fig. [8] As one can see, 
solenoidality of the velocity field is achieved up to the machine accuracy. With 
regard to the modified velocities (the left panel of Fig. [8j), we should empha- 
sized that since the periodic boundary conditions (not homogenous ones) are 
implemented, the modified velocities are not vanished on . In fact, non-zero 
velocities are induced on the T D by the moving obstacle. However, note that 
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Fig. 8. The modified velocities for the moving square cylinder with periodic bound- 
ary conditions. Left: The periodic boundary conditions induced non-zero velocities 
onTfl . Right: However, the velocity vector is solenoidal within the machine accu- 
racies. 

the overall mass balance is satisfied, since the divergences are negligibly small. 
On the other hand, it should be emphasized that satisfying the solenoidality 




(x-x„)/b 

Fig. 9. Errors in immersed boundary conditions of the moving square cylinder due 
to imposing the solenoidality. The positions are presented with respect to the lower 
left corner of the square (xo,yo). Left: Errors in u\ and ui on the left edge of the 
square cylinder. Right: Errors on the upper edge of the square cylinder. 

changes the immersed velocity boundary conditions. Similar problem was ob- 
served in many other immersed boundary methods for the incompressible flow, 
and motivated emergence of some methods like the multi-direct forcing [2"2"f2"5] , 
or the implicit forcing method [TU]|2TI] in the last years. To elucidate the issue, 
the modified velocities on the edge of the square cylinder are compared with 
the desired boundary conditions in Fig. |9j As one can see, both components 
of the velocity (i.e. u\ and ^2) are changed on all the immersed boundaries. In 
fact this is the cost that we have to pay for working with physical (solenoidal) 
velocities. However, a body of evidence show that these discrepancies elimi- 
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Fig. 10. The modified flow field for a moving cube (with velocity u s = (1,0,0)) in 
the quiescent fluid with periodic boundary conditions. Left: Streamlines and the 
contour plots of u\ on half of the solution domain. Right: Iso-surfaces of 10~ 17 level 
of the divergences of the flow. As one can see, the flow is divergence-free within the 
machine accuracy. 



nates by developing the solution in the next times [T5|22f23] ; and furthermore, 
there are some remedies for overcoming this difficulty [T51l23f20] . 



4-3 Extension to the three-dimensional problems 



Although merely the two-dimensional problems are considered up to now, the 
method may be extended easily to the three-dimensional flow. However, in 
these problems all three components of the vorticity vector present, and con- 
sequently we are faced with three Poisson's problems. Therefore, the method 
seems to be not so efficient (in comparison to e.g. some primitive variables 
formulations of the NSE). Nevertheless, since the boundary conditions, and 
the right hand sides of three Poisson's problems ^ are un-coupled, the prob- 
lems are completely parallelizable. Therefore, the method has the potential of 
being competing, at least in principle. 

In the present work, our numerical experiment is restricted to the peri- 
odic boundary conditions (for which we can use the efficient Fourier-spectral 
solvers). The problem set up (which is an extension to the last moving square 



cylinder problem), is shown in Fig. 10 A moving cube Q s = {\\Xi — 71* 00 < 



10/12871"} with velocity u s = (1,0,0) is placed in a regular solution domain 
D = {\\xi — 7r||oo < 27r}, occupied by a quiescent fluid. The extended vorticity 
vector is obtained and the modified velocities are re-calculated on a 128 3 uni- 
form grid, using the Fourier-spectral method. The results are shown in Fig. 



10 As one can see, the cube velocity is induced to the whole of the solution 



domain, and the modified velocity vector is perfectly solenoidal. 
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4-4 Conclusions 



A method was suggested for construction of (nearly) solenoidal velocity vec- 
tors which satisfy (approximately) the desired immersed velocity boundary 
conditions. The method consists of re-calculation of the velocities from an 
extended vorticity field, obtained from an extended initial non-solenoidal ve- 
locity vector. It was shown that the solenoidality of the modified velocities 
are depended on the regular boundary conditions, and two different bound- 
ary conditions were considered. For the homogeneous boundary conditions, 
the method improves the solenoidality, and for the periodic boundary condi- 
tions (that is, the no-boundary problems), the solenoidality satisfies, up to 
the machine accuracy. The modified solenoidal velocities are no-longer satisfy 
the desired immersed velocity boundary conditions, however (as it is shown in 
other references), the discrepancies are eliminated by developing the solution. 
The method was applied to the two- and three-dimensional problems, and 
effects of the method on the velocity fields were discussed. 
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